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Abstract 



We compare results of four public supersymmetric (SUSY) spectrum codes, 
ISA JET 7.71, SOFTSUSY 1 . 9, SPHENO 2 . 2 . 2 and SUSPECT 2 . 3 to estimate the present- 
day uncertainty in the calculation of the relic density of dark matter in mSUGRA 
models. We find that even for mass differences of about 1% the spread in the ob- 
. tained relic densities can be 10%. In difficult regions of the parameter space, such as 

large tan/3 or large mo, discrepancies in the relic density are much larger. We also 
find important differences in the stau co-annihilation region. We show the impact of 
. these uncertainties on the bounds from WMAP for several scenarios, concentrating 

on the regions of parameter space most relevant for collider phenomenology. We 
^Cj ' also discuss the case of Ao 7^ and the stop co-annihilation region. Moreover, we 

present a web application for the online comparison of the spectrum codes. 



1 Introduction 

> ■ 

Since the extremely precise measurement of the cosmic microwave background by the 
WMAP experiment [1, 2], cosmology has been used to severely constrain models with 
cold dark matter candidates. The prime example are supersymmetric models with R- 
parity conservation where the neutralino LSP (lightest supersymmetric particle) is the 
cold dark matter (see Ref. [3] for a review of SUSY cosmology). Requiring that the model 
provide the right amount of cold dark matter 

0.0945 < VLcDuh 2 < 0.1287 (1) 

at 2ct puts strong constraints on the parameter space of the model, in particular in the 
mSUGRA scenario [4, 5, 6, 7, 8, 9, 10, 11]. Effectively, the relic density of dark matter 
imposes some very specific relations among the parameters of the model. Naturally, the 
question arises how precisely Q 1 is calculated in a supersymmetric model. We there- 
fore revisit the constraints from WMAP in the mSUGRA scenario taking into account 
uncertainties originating from the computation of the SUSY spectrum. In the standard 
approach, the relic density is Q oc 1/ (o~v) , where (o~v) is the thermally averaged cross 
section times the relative velocity of the LSP pair. This thermally averaged effective an- 
nihilation cross section includes a sum over all annihilation channels for the LSP as well as 
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co-annihilation channels involving sparticles that are close in mass to the LSP. The relic 
density then depends on all the parameters of the MSSM (i.e. masses and couplings) that 
enter the different annihilation/co-annihilation channels. To calculate the relevant cross 
sections within the context of a model defined at a high scale, say the GUT scale, one first 
needs to solve the renormalization group equations to obtain the MSSM parameters at 
the SUSY scale. Second, higher-order corrections to the masses and couplings need to be 
calculated. Many public or private spectrum calculators perfom this task. The results are 
then used to calculate in an improved tree-level approximation the effective annihilation 
cross-section of neutralinos and the relic density of dark matter. This kind of top-down 
approach is also the typical method to test high-scale models at the LHC [12]. To address 
the issue of the precision of the relic density computation in mSUGRA, in this note we 
compare the results of four public spectrum codes, ISAJET7.71 [13], SOFTSUSY 1 . 9 [14], 
SPHENO 2.2.2 [15] and SUSPECT 2. 3 [16], linking them to micrOMEGAs 1 . 3 . 2 [17] to com- 
pute Q. Since three of these codes, ISAJET7.71, SOFTSUSY 1.9 and SPHENO 2.2.2, are 
of a comparable level as what concerns radiative corrections, the differences in their re- 
sults seem to be a good estimate of the present uncertainties due to higher-order loop 
effects. We also include SUSPECT 2 . 3 in the discussion because it is a widely used pro- 
gram. However, since in contrast to the other three codes, SUSPECT 2 . 3 has only 1-loop 
renormalization group (RG) running for the squark and slepton mass parameters we do 
not use it for the estimate of uncertainties. Within a given program, one can also es- 
timate the theoretical uncertainty by, for example, varying the scale Msusy at which 
the electroweak-symmetry breaking conditions are imposed and the sparticle masses are 
calculated. This was discussed in Ref. [18] and uncertainties on the relic density up to 
20% were found. 

The MSSM parameters that enter the effective annihilation cross section for the LSP 
include all the ones contributing to the annihilation and co-annihilation processes. The 
relic density can then depend on a large number of parameters. However, because one 
needs, at least within the context of SUGRA models, very specific mechanisms to satisfy 
the tight upper bound of WMAP, only a few parameters are critical within each scenario 
[18]. Any shift in one of the critical variable can have a large impact on the value of the 
relic density. Within mSUGRA, the preferred scenarios are the f co-annihilation, the rapid 
Higgs annihilation and the higgsino-LSP scenarios. The main channels are annihilation of 
neutralinos into fermion pairs via s-channel Z or Higgs exchange, or via t-channel sfermion 
exchange, as well as co-annihilation with sleptons. For example, in the co-annihilation 
region, co-annihilation processes are suppressed by a factor exp~~ AM ^ Tf where AM is the 
mass difference with the LSP and Tj m m^o/25 is the decoupling temperature. Then it is 
the mass difference between the NLSP and LSP that introduces the largest uncertainty in 
the prediction of the relic density. In Ref. [18] it was shown that a 1 GeV correction to the 
mass difference could lead to 10% correction on the relic density. In [19] it has been pointed 
out that typical differences in the masses obtained by the spectrum calculator codes are 
of (9(1%), large enough for the computational uncertainty to exceed the experimental one 
of WMAP. In other scenarios, the ones where annihilation proceeds through s-channel 
Z or Higgs exchange, the important parameters are the coupling of neutralinos to the Z 
or Higgs and the mass of the LSP in relation with the mass of the resonance, in general 
the mass of the pseudoscalar. These processes are often relevant in the same "tricky" 
region of parameter space where the discrepancies in the predictions of the spectrum 
calculators well exceed the 1% level, Ref. [19], leading to large uncertainties in the relic 
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density prediction. 

The influence of these differences on relic density computations has first been studied 
in [20] for the Les Houches 2003 workshop. Since then, all above mentioned programs 
have undergone major updates; a re-analysis of the existing uncertainties therefore seems 
appropriate. Moreover, the study of [20] concentrated on potentially large differences 
along specific lines in the focus point, large tan/5 and co-annihilation regions. In this 
article, we consider the WMAP allowed parameter region in the m — mi/2 plane, investi- 
gating in particular differences in WMAP constraints which arise from the different SUSY 
spectrum codes. We also address the issue of non-zero A , which for very large A < 
leads to t co-annihilation. 

We first briefly discuss in Section 2 the calculation of the supersymmetric spectrum. 
We then study in Section 3 some specific scenarios: in Section 3.1 we discuss the case of 
moderate parameters (small m , small to medium m^, moderate tan/3), which is most 
promising for collider phenomenology and where the calculations are expected to be quite 
precise. As it turns out there are, however, non-negligible uncertainties already in this 
region. In Section 3.2, we discuss the case of large tan/5, where much larger differences 
are observed. Section 3.3 then deals with the case of large mo and Section 3.4 with the 
case of large m and large tan (3. Here very large uncertainties are found; in particular 
focus point behaviour may or may not occur depending on the program. The influence of 
the Aq parameter is discussed in Section 3.5. In Section 4, we present a web application 
for online spectrum comparisons. Finally, Section 5 contains conclusions and an outlook. 

For the sake of a fair comparison, we use the same Standard Model (SM) input pa- 
rameters in all programs. In particular, we use mb(mb) MS = 4.214 GeV and a s (Mz) MS = 
0.1172 according to ISAJET7.71. Moreover, we use a top pole mass of m t = 175 GeV 
throughout the paper. The parameters of the MSSM are defined following the SUSY Les 
Houches Accord (SLHA) [21]. 

We do not discuss here the impact of different cosmological scenarii. We assume the 
standard cosmological scenario, in particular that at the freeze-out temperature when the 
interaction rate of particles drops below the expansion rate of the universe, the universe 
was radiation dominated. Modifications of the standard picture for the expansion of the 
universe could significantly affect the estimation of the relic density, examples are models 
with a low-reheating temperature [22] or with scalar- field kination [23]. 

2 SUSY spectrum and relic density 

To derive the relic density within a specific SUSY model, mSUGRA for instance, one needs 
to compute the mass spectrum and couplings from high-energy input parameters. We use 
the latest version of the four public codes ISAJET7.71, S0FTSUSY 1 . 9, SPHEN0 2.2.2 
and SUSPECT 2 . 3 for this task and compare their spectra and the resulting neutralino 
relic densities. These codes basically work as follows: after specifying the gauge and 
Yukawa couplings in the DR scheme at the electroweak scale and starting with an initial 
guess of the MSSM parameters, renormalization group (RG) equations are used to run 
the parameters to some high scale Mx- There boundary conditions are imposed on the 
SUSY-breaking parameters, and the couplings and parameters are run down to the SUSY 
mass scale. At that scale radiative electroweak symmetry breaking is checked. The SUSY 
spectrum is calculated and radiative corrections are computed. The process is repeated 
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iteratively until a stable solution is found. The four programs differ, however, in the 
implementation of radiative corrections (a detailed comparison of the codes can be found 
in [19]). For one, ISAJET7.71, SOFTSUSY 1 . 9 2 and SPHENO 2.2.2 apply full 2-loop RG 
running for all SUSY mass parameters, while SUSPECT 2 . 3 calculates gaugino and Higgs 
mass parameters at 2-loops but squark and slepton parameters only at 1-loop. Second, 
ISA JET 7.71 uses step beta functions when passing thresholds in the RG evolution, adding 
additional finite corrections at the end. In contrast to that the other programs compute 
the complete 1-loop threshold corrections at the SUSY mass scale M S usy = y/ m i 1 m i 2 - 
Third, the use of either on-shell or running masses in the loops can significantly influence 
the results even though the difference is formally a higher-order effect. Moreover, different 
approximations are used in some parts of the loop corrections. For example, ISAJET7.71 
and SPHENO 2.2.2 apply the complete 1-loop corrections given in [24] for the neutralino 
and stau masses, while SOFTSUSY 1 . 9 and SUSPECT 2 . 3 use the approximate expressions 
of [24] for neutralinos and do not include the self-energies for the staus. The calculation of 
the light Higgs mass has recently been standardized between SOFTSUSY 1 . 9, SPHENO 2.2.2 
and SUSPECT 2.3 to full 1-loop plus leading 2-loop corrections, see [25]. ISA JET 7. 71 on 
the other hand uses an 1-loop effective potential, which typically leads to about 2-3 GeV 
higher h° masses compared to the other programs. Notice, however, that this lies within 
the present ~2-3 GeV theoretical uncertainty in rrih- Moreover, as we will see, the exact 
value of rrih is only important in a narrow strip in the large tuq region. All considered, 
we take ISAJET 7 . 71, SOFTSUSY 1 . 9 and SPHENO 2.2.2 as being of a comparable level of 
sophistication as concerns the SUSY and heavy Higgs masses. Two-loop as opposed to 
one-loop scalar running, as in SUSPECT 2 . 3, can however have an important influence on 
the relic density through differences in the sfermion masses. 

The nature of the LSP, which is a linear combination of the bino B, wino W and the 
two higgsino states H\ t 2, is a crucial parameter in the evaluation of the relic density; 

X\ = N U B + N 12 W + iV 13 #i + N U H 2 (2) 

where N is the neutralino mixing matrix. The LSP-higgsino fraction is given as 

f H = N 2 3 + N 2 U (3) 

and is large when the higgsino mass parameter /i < .\/| . where M\ and M 2 are the 
U(l) and SU(2) gaugino masses. The LSP coupling to the pseudoscalar, g^x?^' depends 
on the same elements of the neutralino mixing matrix: 

g mA oc N 2 3 - N 2 U . (4) 

The value of fi depends sensitively, in certain regions of parameter space, on the SM 
input parameters, in particular the top quark mass and its relation with the top Yukawa 
couplings. At large m , the top Yukawa coupling has a strong influence on m 2 H2 ; as a 
result the parameter becomes very sensitive to h t = \/2rht/v 2 , where rh t is the running 
t-quark mass and v 2 is the vev of the second Higgs doublet: 

_ im 2 Hi -m 2 H2 tan 2 g) 1 
~ tan 2 /3-l 2 Z ' [b) 



2 Here note that the default option in SOFTSUSY 1 . 9 is 1-loop running of the squark and slepton mass 
parameters; 2-loop running of these parameters has to be switched on by hand. In the following, we 
always take SOFTSUSY 1 . 9 with full 2-loop RGE running. 



.1 



Here rn 2 H . = m 2 H . — U/vi, % = 1,2, with i{ the tadpole contributions. See [19] for more 
detail. For the intermediate to large values of tan (3 that we will consider, the m 2 H term 
dominates in the extraction of //. Differences in the /i parameter will affect the neutralino 
couplings, in particular the coupling to the pseudoscalar A. In the mass spectrum these 
differences most obviously show up as differences in the mass of the neutralino that is 
dominantly higgsino, usually %\- The programs under consideration all apply the 1-loop 
corrections of [24] plus the 2- loop QCD corrections of [26]. Nevertheless the differences 
in h t are large enough to lead to huge discrepancies in \i at large m . 

The mass of the pseudoscalar, tua, is another important parameter in the computation 
of the relic density. This mass also depends sensitively on the SM input parameters, in 
particular the bottom quark mass and its translation to the bottom Yukawa coupling. 
The bottom Yukawa coupling which is large at high tan j3 impacts the Higgs sector since 
rrVjj is driven by hf, = y/2rhf,/vi, where rh^ is the running b-quark mass and v\ is the vev 
of the first Higgs doublet. The physical pseudoscalar mass directly depends on tt? h : 

The four spectrum codes all apply the corrections of [24, 26], resumming the 1-loop SUSY 
corrections according to [27]. This brings in general good agreement on tua', however as 
we will see the remaining differences can still lead to sizable discrepancies in Q in parts 
of the parameter space. 



3 Results 

3.1 Small mo, small to medium mi/ 2 , moderate tan (3 

We start out with an easy, collider-friendly scenario of small mo, small to medium m\i% 
and moderate tan (3. Such a scenario has gluinos and squarks with masses up to 1 TeV 
which cascade-decay into neutralinos and sleptons. It can hence provide the favourite 
LHC signature of jets plus same-flavor opposite-sign leptons. It also has gauginos and 
sleptons within the kinematical reach of a future e + e~ linear collider (ILC) and is thus 
very well suited for both LHC and ILC studies. 

In the region considered in this section, as in most of the mSUGRA parameter space, 
the LSP is nearly a pure bino. As such it couples preferably to right-chiral sfermions with 
a coupling proportional to the hypercharge. The main annihilation channel for the LSP 
is then into lepton pairs via t-channel exchange of right-chiral sleptons. This process is 
efficient enough to meet the WMAP upper limit only in the low vtlq-vhx^ corner of the 
parameter space, the so-called bulk region. Indeed, for a pure bino LSP the relic density 
is approximately Q oc m~ /tt&o, implying that both the Ir and the Xi must be light. 
Since sleptons must be beyond the reach of LEP2, the upper limit from WMAP is only 
satisfied in a very small region below m\ii ~ 240 GeV, see Fig. 1. The bulk region is, 
however, associated with a light Higgs below the LEP2 limit 3 . Light neutralinos can also 
annihilate efficiently into fermion pairs near a Z or Higgs resonance. This corresponds to 
the near vertical WMAP line in Fig. 1. This possibility is however by large excluded by 

3 Note that an increase in the top-quark mass loosens the LEP2 constraint from the light Higgs. 
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Figure 1: Comparison of results in the m -mi/2 plane, for A = 0, tan/5 = 10, fi > 0, and 
m t = 175 GeV. The red (dark) and orange (light) full lines show the variation of the 2a 
upper limit Q < 0.1287 when micrOMEGAs 1 . 3 . 2 is linked to ISAJET7.71, S0FTSUSY 1 . 9 
or SPHEN0 2.2.2. The orange line basically comes from S0FTSUSY 1 . 9 while the red one 
comes from ISAJET7.71. In addition, the upper bound from SPHEN0 2.2.2 is shown as 
green dotted line, and that of SUSPECT 2.3 as blue dashed line. The light, medium and 
dark gray shaded areas show the regions where the relative differences in Q, 5Q of eq. (7), 
are 4-10%, 10-30% and >30%, respectively. Also shown are contours of minimal (full 
black lines) and maximal (dashed black lines) h° masses as obtained by the spectrum 
codes. The yellow region on the left is excluded by LEP2 constraints; in the yellow 
triangle in the bottom right corner < m^o in ISAJET7.71. The yellow lines show the 
boundaries of the excluded region in the other codes. 



the LEP direct limits [28] on chargino pairs, which in effect translate into a lower limit 
on the LSP mass in mSUGRA. 

Agreement with WMAP is recovered for heavier neutralinos (m^ •> 240 GeV) with 
the additional contributions from co-annihilation channels, the so-called co-annihilation 
region. For co-annihilation to be effective, the mass difference between the slepton NLSP 
and the Xi LSP must be rather small (less than ~10 GeV). Such degenerate sleptons/neu- 
tralinos are found in the low mo region of mSUGRA. The f\ is the lightest slepton due both 
to the effect of the r- Yukawa coupling in the RGE running of rrif R as well as to the mixing 
between tl and tr which lowers the mass of the fi to m fl < rrif R . In fact, co- annihilation 
processes with fi dominate over most of the allowed region in Fig. 1. In the co-annihilation 
region, it was shown in Ref. [18] that the relic density is extremely sensitive to the mass 
difference between f\ and x?> AM^^i)- Typically a shift in AM(xSfi) pa 1 GeV induces 
AQ pa 10%. Previous comparisons between the public SUSY spectrum codes [19] have 
shown that the predicted masses often differ by more than ±1%, inducing discrepancies 
in AM(xiTi) above ±1 GeV and hence large uncertainties in the relic density. 

These expectations are corroborated by a scan in the mo-mi/2 plane comparing the 
predictions of the four spectrum codes. Figure 1 shows results for A = 0, tan/5 = 10, 
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[i > and rat = 175 GeV. The red and orange lines show the variation of the 2a upper 
limit Q < 0.1287 when micrOMEGAs 1.3.2 is linked to ISAJET7.71, SOFTSUSY 1 . 9 and 
SPHENO 2.2.2. In addition, the light, medium and dark gray shaded areas show the 
regions where the relative differences in Q, 

= iS^max ^min) /^mearu (7) 

are 4-10%, 10-30% and >30%, respectively. Here Q max and fi min are the maximal and 
minimal values and VL mean the arithmetic mean of the Q values obtained from ISA JET 7 . 71, 
SOFTSUSY 1 . 9 and SPHENO 2 . 2 . 2 at a specific parameter point. We do not include SUSPECT 2 . 3 
in the calculation of 5q because it has only 1-loop scalar running. A 5q of 30% corresponds 
to the present precision of WMAP, while the PLANCK experiment [29] is expected to 
reach a precision of 4%, corresponding to the white area in Fig. 1. Also shown are con- 
tours of the minimal and maximal h° masses as obtained by the four spectrum codes. As 
a general rule, the m™ ax lines come from ISAJET7.71, while the m™ m lines come from 
the other programs. Note that the bulk region is practically excluded. 

The red (maximal Q) and orange (minimal Q) lines in Fig. 1 come from ISAJET7.71 
and SOFTSUSY 1 . 9, respectively. The values obtained from SPHENO 2.2.2, shown as dot- 
ted green line, lie in between these curves. In the co-annihilation region the results of 
SUSPECT 2 . 3, shown as dashed blue line, fall within the red and orange lines. While the 
differences in the WMAP bounds in Fig. 1 do not look dramatic, it becomes clear from 
the grey shaded areas that the relative differences in Q are quite large in the allowed 
parameter space, that is in the co-annihilation region, where the precise mass differences, 
in particular between fi and x?; are important. 

ISAJET 7 . 71, SOFTSUSY 1 . 9 and SPHENO 2.2.2 typically agree on the t\ mass to ~1%. 
The difference mainly comes from SOFTSUSY 1.9, which neglects the tau Yukawa coupling 
h T and the f self-energy correction and hence gets a slightly smaller m^. SOFTSUSY 1 .9, 
SPHENO 2.2.2 and SUSPECT 2 . 3 agree very well on the xi mass, while ISAJET 7.71 finds a 
m^o smaller by about 2%. As a consequence, both SOFTSUSY 1 . 9 and SPHENO 2.2.2 tend 
to give a smaller fi—Xi mass difference than the other two programs, and hence a smaller 
Q in the co-annihilation region. As an example, Table 1 lists the relevant masses together 
with fl for m = 70 GeV, mi/2 = 350 GeV, A = 0, tan/5 = 10 and /i > 0. Table 2 gives 
the according relative contributions to Q^ 1 for this point. Note here the contribution of the 
fi and eR, JIr co-annihilation channels. Clearly our expectations that the mass difference 
is the most important parameter are confirmed. The 2 GeV decrease in AM(x5fi) when 
going from SPHENO 2.2.2 to SOFTSUSY 1.9 roughly corresponds to a decrease of 0(20%) 
in Q as expected. As a result of the mass spectrum, one finds a larger contribution from 
the co-annihilation channels for SOFTSUSY 1 . 9 where it amounts to almost 80% of the 
effective annihilation cross section as compared to the other codes where co-annihilation 
channels contribute ~50-70%. ISAJET 7.71, which agrees well with SPHENO 2 . 2 . 2 on the 
fi mass but finds a smaller m^o, has the largest Xi^i mass difference and a ~50% higher 
Q as compared to SPHENO 2.2.2. For similar AM(%?fi), ISAJET 7. 71 predicts a slightly 
lower value for the relic density as compare to other codes because of a lower LSP mass. 
SUSPECT 2 . 3 on the other hand agrees well with SOFTSUSY 1 . 9/SPHENO 2 . 2 . 2 on the LSP 
mass, but due to the missing 2-loop effects in the running of the slepton masses it gets 
a heavier f\ (and e#) and hence a larger Q. We have checked that when using only 1- 
loop RGEs for the slepton mass parameters in SOFTSUSY 1.9, it reproduces the results of 
SUSPECT 2 . 3. 
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ISAJET7.71 


SOFTSUSY 1 . 9 


SPHENO 2.2.2 


SUSPECT 2 . 3 


X? 

f\ 

en 
h° 


136.7 
147.7 
155.7 
115.8 


140.0 
145.7 
153.8 
113.1 


139.5 
147.1 
155.4 
113.4 


140.0 
149.7 
157.6 
113.3 


mfj — m^o 


11.0 


5.7 


7.6 


9.7 


n 


0.136 


0.069 


0.092 


0.120 



Table 1: Relevant masses, the x?-fi mass difference (in GeV) and the resulting Q for 
m = 70 GeV, mi/2 = 350 GeV, A = 0, tan/3 = 10 and \i > 0. The higgsino fraction of 
x\ is 1.4-1.5% in all cases. 



channel 



ISAJET7.71 SOFTSUSY 1.9 SPHENO 2.2.2 SUSPECT 2. 3 



Xi T i - 
hn - 
nfi - 



■> ee 
■> rr 

+ 7 /Ze 

► 7/Zr 
rr 

77, 7 Z 

► er 



28% 
16% 
8% 
30% 
5% 
2% 
2% 



10% 

6% 
8% 
39% 
17% 
7% 
6% 



16% 
9% 
8% 
38% 
11% 
6% 



22% 
13% 
10% 
34% 
7% 
3% 
2% 



Table 2: Relative contributions to Q 1 for m = 70 GeV, m^j = 350 GeV, Aq = 0, 
tan j3 = 10 and \i > 0; with e = e, n and e# = e^, 

Some more comments are in order. First, a non-zero value of Aq shifts the contours 
of constant Higgs masses and moves the position of the stau co-annihilation strips as well 
as of the excluded regions in Fig. 1; it does however not change the picture qualitatively, 
provided Aq is not so large as to make t\ the (N)NLSP. The case of a non-zero Aq will 
be discussed in detail in Section 3.5. Second, for the reference point SPSla' of the SPA 
project [30], (m = 70 GeV, m 1/2 = 250 GeV, Aq = -300, tan/3 = 10, fi > and 
m t = 178 GeV), ISAJET 7 .71, SOFTSUSY 1 . 9 and SPHENO 2.2.2 give values of fi = 0.126, 
0.103 and 0.114, respectively. SOFTSUSY 1.9 with 1-loop scalar running gives Q = 0.125, 
and SUSPECT 2 . 3 SI = 0.126. All values lie within the WMAP allowed range of Eq. (1) at 
this point, with the spread of Sfl ~ 20% again being mainly due to AM(x5^i)- 

3.2 Large tan (3 

We next consider large values for tan/3; we stay however within collider-friendly scenarios 
with small rriQ and small to medium mu2- At large values of tan /3, the enhanced couplings 
of the heavy Higgses to bb and rr lead to an enhancement of neutralino annihilation 
channels through XiXi A — > bb, tt. Because of the Majorana nature of the 

LSP the main contribution is the pseudoscalar exchange, the CP-even state being P- 
wave suppressed. Even though the LSP is mostly bino, its small higgsino component is 
sufficient to make annihilation through the pseudoscalar and the Goldstone component of 
Z exchange dominant. These contributions are added to the contributions from t-channel 
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sfermion exchange or from co-annihilation with staus that were already present at lower 
values of tan (3. 

At low mi/2, the annihilation into bb typically constitutes more than 80% of the effec- 
tive annihilation cross section. For a fixed value of mi/2, hence of neutralino mass, the 
relic density decreases with mo since both the sfermion masses as well as the pseudoscalar 
mass decrease, making for more efficient annihilation. Because of the enhanced contribu- 
tion of the pseudoscalar exchange, a much larger region of parameter space in the bulk 
is compatible with the WMAP upper bound as compared to intermediate tan (3 values, 
see Fig. 2. Nevertheless as one moves towards larger values of m\/2 and a heavier LSP, 
one must again appeal to co-annihilation to retain consistency with WMAP, leading to a 
mixed region with both co-annihilation and pseudoscalar exchange. The co-annihilation 
occurs exclusively with f±, which is much lighter than the other sleptons at large tan/3. 
Note that for the range of m\/2 which we are considering, we are never near the heavy 
Higgs pole. 

The relic density is again sensitive to AM(x\f\) for the co-annihilation processes. 
Sensitivity to AM(xJil) = m^- 2m^o as well as to the X1X1A coupling are expected for the 
Higgs contribution, see [18]. As already mentioned, the bottom and tau Yukawa couplings 
play an important role in radiative corrections to the sparticle and Higgs masses at large 
tan /3, leading to larger differences in the spectra. Consequently in the computation of 
the relic density we also observe larger discrepancies between the four codes. 

Figure 2 compares the results of the various codes in the m -mi/ 2 plane analogous to 
Fig. 1 but for tan/3 = 40 (left) and tan/3 = 50 (right). The other parameters are Aq = 0, 
/i > and m t = 175 GeV as before. At tan/3 = 40, the WMAP exclusion curves seem to 
agree quite well. Small differences (AQ < 5%) are observed over much of the plane, but 
these increase rapidly to 10-30% and more as one moves into the WMAP allowed region. 
Near the stau-LSP border, differences in the predictions of the spectrum calculators for 
the f\ masses, and hence for — m^o, explain this discrepancy, just as was the case for 
tan/3 = 10. Large differences are also observed for low mi/2 in the region near the band 
excluded by LEP limits. These discrepancies are due to differences in m^o. Specifically 
some codes allow a significant annihilation rate through the light Higgs exchange in a 
region that is allowed by the LEP limit on charginos. Here again the low mi/2 (bulk) 
region is not compatible with the lower limit on the Higgs mass. 

For further illustration, we pick a parameter point from Fig. 2a, (mo, m^) = (194, 300) 
GeV at tan/3 = 40. Details on the spectrum relevant for the relic density calculation and 
the list of important channels for all four codes are presented in Table 3. For this pa- 
rameter choice we are in a mixed region where both co-annihilation and Higgs exchange 
processses are important. All codes agree quite well on the values of m^o, m^ and con- 
sequently rriA — 2m^o with maximal variation on the latter of about 5%. The variation 
in the /i parameter is below 3%. The variation of the NLSP-LSP mass difference is 
2% within ISAJET7.71, S0FTSUSY 1 . 9 and SPHEN0 2.2.2, but 40% if we also include 
SUSPECT 2 . 3. The difference in the n mass between S0FTSUSY 1 . 9 and SPHEN0 2 . 2 . 2 can 
again be explained by the missing h T and f self-energy corrections in the former program, 
which is roughly a 1% effect. It is interesting to note that this also influences rriA at the 
level of few per-mille. All considered, the uncertainties in the Higgs annihilation and the 
stau co-annhilation channels are of similar importance in ISAJET 7 .71, S0FTSUSY 1 . 9 and 
SPHEN0 2.2.2, leading to a spread in Q of 25% for the parameter point of Table 3. If we 
interpret this as Q = 0.107 ± 0.013, then SUSPECT 2.3 deviates by 2.7<r due to its larger 
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Figure 2: Comparison of results analogous to Fig. 1 but for tan (3 = 40 (left) and tan /? = 
50 (right); Aq = 0, /i > 0, and m t = 175 GeV. The red and orange lines show again the 
variation of the bound Q < 0.1287 due to differences in the spectra from ISAJET7.71, 
S0FTSUSY 1 . 9 and SPHEN0 2.2.2. In the right panel, the blue dashed line shows in addition 
how the upper curve would move when including SUSPECT 2 . 3. 
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Figure 3: WMAP allowed regions of 0.0945 < Q < 0.1287 for tan/3 = 50, A = 0, 
fj, > 0, m t = 175 GeV; left: ISAJET7.71 and S0FTSUSY 1 . 9, right: SPHEN0 2.2.2 and 
SUSPECT 2 . 3. 
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NLSP-LSP mass difference. 

Increasing further tan/3 means that the Abb coupling is further enhanced and pseu- 
doscalar exchange dominates over most of the probed parameter space. Figure 2b com- 
pares the various codes for tan/3 = 50. In this case, large discrepancies are found in 
the relic density and this over most of the parameter space. In particular, the minimal 
and maximal upper boundaries from WMAP displayed as red and orange curves differ 
significantly. As a result of the pseudoscalar exchange contribution, the bulk region is 
much larger compared to small tan/3. We however stress that this bulk region is not 
one of t-channel sfermion exchange but rather one of heavy Higgs annihilation. Stau co- 
annihilation also plays a role near the stau-NLSP boundary; however this typically drives 
the relic density below the WMAP range. Large tan/3 also means larger discrepancies 
between the predictions of the spectrum calculators especially for the pseudoscalar mass. 
The large discrepancies in Q over most of the parameter space are due to the differences in 
rriA — 2m^o. Typically, at such large tan/3 the mass of the pseudoscalar in SOFTSUSY 1 .9 
and SUSPECT 2 . 3 is lighter than that predicted by the other two codes. The only region 
where 5Q < 10% lies at small mi/2 where the pseudoscalar exchange diagram is less 
important and better agreement for the masses is found. 

Figure 3 explicitly compares the 2a WMAP allowed regions of the four programs 
at tan (3 = 50. The difference in the prediction of the pseudoscalar mass also explains 
why the band near the stau co-annihilation region is much narrower for ISA JET 7.71 and 
SPHENO 2.2.2. Here one is sitting too far from the Higgs resonance to get a significant 
contribution to the annihilation cross- sect ion, the only remaining WMAP allowed region 
being the narrow stau co-annihilation strip. For further illustration, we pick a parameter 
point from Fig. 3, (m , m^) = (350, 350) GeV at tan/3 = 50. The results for this point 
are listed in Table 4. Here the mass difference between the NLSP and the LSP is much too 
large to get a significant contribution from co-annihilation processes. The main channels 
are annihilation of neutralinos into fermion pairs via pseudoscalar exchange. Although 
one is far from the Higgs resonance, this process is efficient enough due to the enhanced 
coupling of the Higggs to bb and rr. For this point one gets as usual rather good agreement 
among all codes in the x? masses and in the higgsino fraction. The pseudoscalar masses 
also agree within 1-2%; for the resonance parameter mA—^m^o the differences are however 
around 10%. In [18] it was shown that in this region a 4% shift in tua — 2m^o leads to 
a 10% change in Q. The discrepancies in the mass difference found in Table 4 explain 
the difference between the value of the relic density in SPHENO 2.2.2, SOFTSUSY 1 . 9 and 
SUSPECT 2.3. In the case of ISA JET 7.71 the decrease in the annihilation cross section due 
to the fact that one is sitting further away from the Higgs resonance is partly compensated 
by a lower value of the /i parameter (to wit the smaller value of m^o) hence a larger 
XiX^A coupling. We have also checked explicitly that by adjusting rriA — 2m^o to the 
SPHENO 2.2.2 value in SOFTSUSY 1 . 9 we recover very good agreement between the two 
programs. 

To put these results in perspective, we also remark that there is a strong mb(mb) 
dependence in the computation of the pseudoscalar Higgs mass as discussed in Section 2. 
This has an impact on the relic density [31, 18]. For example for the parameters of 
Table 4, decreasing m^rrib) to 4.168 GeV (less than a 2% change) makes the result of 
SOFTSUSY 1 . 9 agree perfectly with the ones from ISAJET 7 . 71. Considering that there are 
large theoretical uncertainties in the extraction of m^m^), this source of uncertainty at 
present exceeds the one estimated by taking the difference between codes. 

i l 





ISAJET7.71 


S0FTSUSY 1 . 9 


SPHEN0 2.2.2 


SUSPECT 2 . 3 


Xi 


117 9 
11 1 .Z 


1 1 Q Q 

i iy .y 


1 1 Q 7 

i iy. i 


1 1 Q Q 

i iy.y 


fi 


lO 1 .41 


1 11 9 


1 11 A 
101.41 


1 "37 7 
lO ( . 1 


n 


111;') 
110. o 


119 7 
1 1Z. 1 


1 lo.U 


1 1 9 8 
llZ.o 


A 


OUO.'i 


ouo.z 


1P.P, A 




v° 

A3 


394 Q 


401 4 


40 s ) 3 


40 s ) 3 


777 ~ — 777 „ n 


14.2 


13.3 


11.6 


17.8 


rriA — 2m^o 


129 


123 


127 


125 


XiXi - bb 


40% 


38% 


30% 


49% 


X°iXi -> ee 


12% 


10% 


10% 


14% 


XiXi -> 


17% 


14% 


13% 


19% 


X?fi -> /it 


13% 


16% 


21% 


7% 




12% 


14% 


18% 


7% 


fifi — > /i/i 


1% 


2% 


3% 




fi 


0.120 


0.107 


0.094 


0.142 



Table 3: Masses and mass differences (in GeV), the most important contributions, and 
the resulting fi for m = 194 GeV, m 1/2 = 300 GeV, A = 0, fi > and tan/5 = 40. The 
higgsino fraction of Xi is 1-8% in all cases. 
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Table 4: Same as Table 3 but for rriQ 
fraction of Xi is l- z 
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3.3 Large m , focus point 



Large tuq is a notoriously difficult region which suffers from large uncertainties. The reason 
is the extreme sensitivity of the /i parameter to the top Yukawa coupling alluded to in 
Section 2. We limit our discussion to gaugino and higgsino masses within the reach of LHC 
and ILC and consider values of mo up to 4.5 TeV. Figure 4 shows the allowed regions in 
the rrio-min plane for m = 1-4.5 TeV, tan (3 = 10, A = 0, \i > and m t = 175 GeV. A 
striking discrepancy between the codes is the occurrence or non-occurrence of focus-point 
behaviour and related with this the limit of radiative electroweak symmetry breaking 
(REWSB). In Fig. 4, the four programs agree more or less up to mo ~ 2 TeV. Above 
this value, the results of ISAJET 7 .71 become very different, with REWSB breaking down 
around m ~ 2.7-3 TeV. In SOFTSUSY 1.9 and in SUSPECT 2. 3 this happens only around 
rrio ~ 3.5-4 TeV while in SPHENO 2.2.2 one can go to much higher mo- In fact this 
behaviour is related to small differences in the treatment of the top Yukawa coupling; 
focus point behaviour can be recovered for all codes when one lowers the top-quark mass. 

In the allowed parameter space, the main annihilation channel for neutralinos is into 
fermion pairs. Consistency with WMAP then requires some enhancement factor for the 
annihilation cross section. This is in principle provided by the light Higgs resonance - 
but only in a narrow strip of the parameter space. The relic density hence becomes very 
sensitive to the m^ — 2m^o mass difference (as compared to the decoupling temperature 
of the neutralinos, Tj pa m^o/25). The width of the h° is not an important parameter 
because it is much smaller then Tf. 

When the Xi mass is slightly below half the h° mass, most of the x?' s annihilate effi- 
ciently through XiXi ~^ ^° ~^ bb. This requires a very small mi/2, roughly mi/2 < 150 GeV 
as can be seen in Fig. 4. On the other hand, the LEP bound of m-± > 103 GeV requires 
mu2 > 130-140 GeV. In Fig. 4 the bands that are within the 2a WMAP range corre- 
spond to rrih — 2m^o either of a few hundred MeV or around 10 GeV (15 GeV in case of 
ISAJET 7. 71). In between these values, the Higgs annihilation mechanism is too efficient, 
resulting in Q < 0.0945. Table 5 gives examples for mo = 2 TeV and mo = 3.8 TeV. For 
m = 2 TeV, SOFTSUSY 1 . 9, SPHENO 2.2.2 and SUSPECT 2 . 3 predict similar masses and 
LSP higgsino fractions. As expected, the relic density decreases as one moves slightly away 
from the pole. For ISAJET 7. 71, predictions for the relic density for a given — 2m^o 
are typically lower than for the other codes, since two other effects enhance the annihi- 
lation cross-section: a larger LSP higgsino fraction and the fact that with a lighter LSP 
one benefits from the Z-exchange contribution. For m = 3.8 TeV, only SOFTSUSY 1.9, 
SPHENO 2.2.2 and SUSPECT 2 . 3 find viable RGE solutions. There are now large ~ 60% 
discrepancies in the \x parameter also among these codes. This is reflected in quite dif- 
ferent higgsino fractions, and in turn in (9(100%) differences in the values of Q. As a 
side remark we note that the large uncertainties in the \x parameter also lead to signifi- 
cant discrepancies in the X3 4 an d X12 masses, which can considerably impact the collider 
phenomenology of a particular mSUGRA point. 

Another comment is in order. Within any of the spectrum codes a change in m t of the 
order of what will be measured at LHC (Am* ~ 1 GeV) induces large changes in the value 
of fi and hence in the LSP mass, its higgsino fraction, and the relic density. The latter 
can vary by over an order of magnitude within a given code. This is due to the extreme 
sensitivity of the running of m 2 H2 to the top Yukawa coupling as explained in Section 2, 
c.f. eq. (5). Small changes in the input value of m t can therefore bring approximate 
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Figure 4: WMAP allowed regions (blue) in the m -mi/2 plane for large m ; tan/? = 10, 
A = 0, fi > 0, m t = 175 GeV. In the dark blue bands 0.0945 < Q < 0.1287, while in the 
light blue bands Q < 0.0945. In the gray areas there is no radiative EWSB; the yellow 
regions are excluded by the LEP bound m-± > 103 GeV. 
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Table 5: Relevant masses and mass differences (in GeV), the higgsino fraction of the 
LSP, the most important contributions and the resulting Q for mo = 2 and 3.8 TeV, 
mi/2 = 144 GeV, and tan/3 = 10 (A = 0, // > 0, m t = 175 GeV). 

agreement between the different codes. We emphazise however that this only reflects the 
large theoretical uncertainty in this regime. 

3.4 Large mo, large tan (3 

As we increase tan f3 it becomes increasingly easier to reach the focus point region. There 
is also a strong dependence on the value of the top-quark mass, and typically ISAJET 7.71 
can find a focus point behaviour with significantly heavier m t than the other codes [20]. 
We consider in more details the case tan/5 = 50 and m t = 175 GeV. A value for the relic 
density in agreement with WMAP requires Mi < \i < M2 so that the LSP is a mixed 
bino-higgsino state. As one moves very close to the electroweak symmetry breaking border 
and fi drops even below Mi, the higgsino fraction increases rapidly; the relic density 
drops below the WMAP range. In what follows we concentrate again on collider-friendly 
scenarios with not so heavy neutralinos and charginos. 

At large tan j3 and large tuq, the main neutralino annihilation channels are into fermion 
pairs or into pairs of gauge or Higgs bosons. Fermion pair production proceeds through 
s-channel exchange of Higgs or Z (the Goldstone component) and is proportionnal to 
the fermion mass. Annihilation into tt is therefore favoured as soon as it becomes kine- 
matically accessible. If not, W-pair production is the dominant channel, proceeding via 
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Table 6: Relevant masses (in GeV), the higgsino fraction of the LSP, the most important 
contributions and the resulting Q for m = 3450 GeV, m\n = 350 GeV, tan/3 = 50, 
A = 0, (i > 0. 

t-channel exchange of charginos. Neutralino/chargino co-annihilation channels can be im- 
portant as well, but typically they are so efficient that they lead to Q below the WMAP 
range. The \x parameter, which determines the neutralino and chargino masses as well as 
the XiXi^/A coupling, also has a significant influence on the relic density. In [18] it was 
shown that A/x ~ 1-2% could induce shifts of 10% in fl for tan/3 = 50. The dependence 
on m^o or on the pseudoscalar mass is expected to be weaker; corrections of 50% or larger 
are necessary to induce a 10% shift in Q [18]. 

We consider more closely the point mo = 3450 GeV, mi/j = 350 GeV and tan/3 = 50. 
Table 6 displays the results for the spectrum and the most important contributions. For 
this scenario S0FTSUSY 1 . 9 gives a result within the WMAP range. ISAJET 7 . 71, however, 
does not find a solution to the RGEs; the relic density of SPHEN0 2.2.2 and SUSPECT 2 . 3 
is orders of magnitude above the WMAP bound. The reason for the latter is that the 
fi parameter in SUSPECT 2 . 3 and even more in SPHEN0 2.2.2 is much larger and one is 
in a regime of a mostly bino LSP -hence no efficient channel for annihilation is available 
(c.f. the values for m^o and in Table 6). Moreover, S0FTSUSY 1.9 predicts much 

lighter charginos, which makes annihilation into W pairs through chargino exchange more 
efficient. Owing to the huge Xi~A° mass difference, the influence of the pseudoscalar mass 
on the relic density is small, although there is a spread in the prediction of vtla of several 
hundred GeV. This discrepancy in tua becomes, however, very relevant in the Higgs funnel 
region, that is for larger values of mm. 

As was the case in the previous section, within any of the spectrum codes a small 
change in m t induces large changes in the value of \x and hence the relic density, which 
can vary as before by over an order of magnitude within a given code. Using a different 
input value for mt can therefore compensate the large discrepancies observed between 
different codes. For example, a decrease of about 0.5 GeV in m t brings the results of 
SUSPECT 2 . 3 for both the spectrum and the relic density, in good agreement with those of 



SOFTSUSY 1.9. SOFTSUSY 1 . 9's results of Table 6 can also be approximately reproduced 
with ISAJET7.71 using m t = 176.36 GeV. Note however that this amounts to extreme 
fine-tuning. 

3.5 Varying A 

Non-zero values of A can significantly influence the scalar masses as well as the fi pa- 
rameter. Roughly speaking, for A < 4 the ti,&i,fi masses decrease while m A and \i 
increase. For A > 0, the shifts go in the opposite directions. The pseudoscalar mass is 
relevant for annihilation processes at large tan f3 where XiXi A ^ ff. The \x param- 
eter determines the higgsino fraction of the LSP. It also directly influences the mixing in 
the stau sector and therefore the contribution of the f coannihilation processes. With 
our convention, A = leads to A t < at the weak scale; A < increases \A t \ (at 
the weak scale) thus lowering the t\ mass through i) RG running and ii) a larger tir^n 
mixing. Analogous arguments hold for sbottoms and staus, though here the L R mixing 
is dominated by \x tan (3. Also the running of A T is less strong, so that A T usually does not 
change sign with respect to Aq. The masses and trilinear couplings of the third generation 
enter in turn the running of the Higgs mass parameters, the radiative corrections to the 
Higgs pole masses, and the computation of \i. 

The uncertainties in the masses, estimated as the differences between the codes, tend 
to be larger for A ^ as compared to A = 0. Nevertheless the general picture outlined 
in the previous sections holds, as the same mechanisms as for Aq = are at work for 
neutralino (co-) annihilation over most of the parameter space. Only when \A \ becomes 
large enough to make t\ very light, in fact the NLSP or NNLSP, new co- annihilation 
channels appear associated with a new region of parameter space where the relic density 
is consistent with WMAP. 

Let us discuss the cases of moderate and large tan/5 in more detail. For tan/5 = 10, 
a non-zero value of A shifts the contours of constant light Higgs masses (towards lower 
values of min for A < 0) and moves the position of the stau co-annihilation strips as 
well as of the excluded regions (towards higher values of m for A < 0) as compared to 
Fig. 1. It does however not change the picture qualitatively; the WMAP allowed regions 
are a small bulk region with XiX? annihilation and a narrow strip of co-annihilation with 
staus. There is an increase in the differences in m~ — m^o and hence in Q between the 
codes, but the effect is in general not very large. The only new feature appears for values 
of \Aq\ large enough to make ti the (N)NLSP. This case will be illustrated later in this 
section. 

For tan (3 = 50, we observe larger discrepancies between the codes even for moderate 
values of Aq. This is not surprising as the mixing in the f sector depends on /itan/5 
and relatively small shifts in \i can have important effects on the f\ mass. Moreover, the 
pseudoscalar mass is quite sensitive to A . While for A = (and small to medium 
mo) ISAJET7.71, SOFTSUSY 1.9, SPHEN0 2.2.2 and SUSPECT 2.3 typically agree on m A 
to 1-2%, for Aq ^ differences of a few per-cent can show up. Figure 5 compares the 
regions of the mo-wii/2 plane compatible with the upper limit of WMAP analogous to 
Fig. 2b but for A = mm. As can be seen, the WMAP bound is shifted towards higher 
values of m . This is because, as mentionned above, the pseudoscalar mass decreases 

4 Using SLHA conventions, the off-diagonal element of the {up, down}-type sfermion mass matrix is 
m 2 LR = (Af — /i{cot/3,tan/3})m/. 
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with increasing Aq, so annihilation channels through Higgs exchange are favoured. The 
Higgs exchange dominates over most of the region of the plot, however with rather large 
differences in Q. The largest differences are found for m^j > 400 GeV, as was the 
case for Aq = 0. Differences in the pseudoscalar masses increase with increasing Aq and 
mi/2, with SOFTSUSY 1 . 9 and SUSPECT 2 . 3 predicting smaller m A than SPHENO 2.2.2 and 
ISAJET7.71 for A = my 2 > 300 GeV (for smaller values, it is ISAJET7.71, which 
predicts the lightest m^). ISAJET7.71 also predicts a lighter LSP and a lighter fi in 
the co-annihilation range and hence a much lower Q as one moves closer to the fi-LSP 
boundary. SUSPECT 2 . 3 on the other hand predicts larger f\ but smaller pseudoscalar 
masses than the other programs. This leads to a larger value for Q from the SUSPECT 2 . 3 
spectrum for Aq = m 1( / 2 > 400 GeV. In Fig. 5, the dashed blue line shows how the 
exclusion curve correponding to the maximal Q is shifted when including SUSPECT 2 . 3. 
One can conclude that for tan /3 = 50 differences between the codes are large everywhere, 
with 5Q exceeding 30% in a large portion of parameter space. 

We have also studied the case Aq = —mi/2 at tan/3 = 50. Here the pseudoscalar 
mass is larger than in the Aq = case, so a relic density in agreement with WMAP 
requires, especially at large m^, some contribution from co- annihilation processes, in 
particular with f\. Therefore the value of the relic density is once again very sensitive 
to the X?-Ti mass difference, and discrepancies in Q are larger than for A = 0. Since 
for Aq < we encounter instabilities in the scan with SPHENO 2.2.2, we do not show a 
plot but examplify this case in Table 7 for m = 376 GeV, mi/ 2 = — Aq = 400 GeV, 
tan/3 = 50 and fi > 0. Here the discrepancy in the f\ mass reaches about 10%, meaning 
that the Xi~^i mass difference varies by more than 100%, thus inducing huge differences 
in the relic density. The lightest f\ is again obtained with ISAJET7.71. We do however 
find rather good agreement between the codes as concerns the boundary of the WMAP 
region. 

A special case is a very large negative Aq, such that t\ becomes light enough to 
contribute to co-annihilations. This is the case when the t\ is the NNLSP or even the 
NLSP. The relic density is then very sensitive to the mass difference between x\ and 
t\. Since the largest discrepancies between spectrum calculators are usually found for 
the masses of coloured sparticles [19], the predictions for the relic density and for the 
region compatible with WMAP can differ significantly in this case. Figure 6 shows the 
WMAP-allowed strips in the rriQ-mi/2 plane for Aq = —Ami/2, tan/3 = 10 and fi > 0. 
For ISAJET7.71, SOFTSUSY 1.9 and SPHENO 2.2.2, co-annihilation with stops dominates 
when rrii/2 < 350-400 GeV, while for larger m^j one has mostly stau co- annihilation. For 
SUSPECT 2 . 3, ti co-annihilation dominates over the whole allowed region. As expected, the 
allowed bands are very narrow. They correspond to — m^o 20-30 GeV, typically a 
much larger mass difference than for the case of f co-annihilation. This is due to the large 
cross section of Xi^i ~ > tK ^9- Table 8 shows the spectrum as well as the most important 
contributions to Q for one point, m = 161 GeV, = 350 GeV, Aq = —1400 GeV, 
tan P = 10 and ji > 0. As one can see, SOFTSUSY 1 . 9 and SPHENO 2.2.2 agree quite well 
on the ti mass and hence the relic density, with only few per-cent difference between the 
two programs. In comparison, ISAJET7.71 predicts a lighter t\ and thus a much smaller 
relic density at a given parameter point. The difference is at the level of 10% for 
and of (9(100%) for Q. Also the boundaries where ti becomes the LSP are quite different 
between SOFTSUSY 1 . 9/SPHENO 2 . 2 . 2 on the one side and ISAJET 7 . 71 on the other side. 
Part of the discrepancies may come from large logs in the RGEs in ISAJET 7. 71 due to 
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Figure 5: Comparison of results analogous to Fig. 2b (tan/3 = 50) but for A = m^. 
The red and orange lines show the variation of the bound ft < 0.1287 due to differences 
in the spectra. The dashed red and orange lines show the situation when only comparing 
S0FTSUSY 1 . 9 and SPHEN0 2.2.2. The gap between the dashed and the full red lines is 
due to a lighter f\ and hence more f co-annihilation in ISAJET7.71; the gap between the 
dashed and the full orange lines is due to smaller Xi an d ^4° masses in ISA JET 7. 71. The 
dashed blue line shows again how the maximal ft moves when including SUSPECT 2.3. 
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fifi -> WW, ZZ, 77 


15% 
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fl 


0.017 


0.107 


0.081 


0.136 



Table 7: Relevant masses and mass differences (in GeV), the most important contributions 
and the resulting ft for m = 376 GeV, m 1/2 = 400 GeV, A = -400 GeV, tan (3 = 50 
and /i > 0. The higgsino fraction /j/(x?) is 0.8%. 
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Figure 6: WMAP strips from the four public codes for A = — Ami/2, tan/5 = 10, \i > 
and m t = 175 GeV. The yellow region in the bottom right corner is excluded due to a f\ 
LSP. In the yellow bottom left region, S0FTSUSY 1.9 and SPHEN0 2.2.2 have a h LSP; 
the yellow dashed line shows the bound of h LSP in ISAJET7.71. 
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Table 8: Relevant masses and mass differences (in GeV), the most important contributions 
and the resulting ft for m = 161 GeV, m 1/2 = 350 GeV, A = -1400 GeV, tan/3 = 10 
and fi > 0. f H (x°i) ^ 0.4%. 
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the very large mass splitting of the stops. Much larger discrepancies are however found 
when comparing with SUSPECT 2 . 3. Since SUSPECT 2 . 3 does not have the 2-loop RGEs 
for the squark parameters, including A t , it predicts a much lighter t\ than the other three 
programs. For the point of Table 8, the difference in mi is about 60 GeV, or 35%, making 
t\ the LSP in the SUSPECT 2 . 3 spectrum. The large discrepancy between SUSPECT 2 . 3 and 
the other programs can be seen clearly in Fig. 6. Again, the results of SUSPECT 2. 3 are 
reproduced by using only 1-loop RGEs for squark and slepton parameters in SOFTSUSY 1 . 9. 
Last but not least notice also that for the parameters of Table 8, in SOFTSUSY 1.9 and 
SPHENO 2.2.2 even though f\ is the NLSP, coannihilation channels with t\ dominate. 

In summary, at very large Aq < one can get phenomenologically very different 
scenarii for the same mSUGRA point; it is clear that including the full two-loop RG 
running plus a careful treatement of threshold corrections is important for a reliable 
prediction of the relic density. 

In this context it is also interesting to compare with the results of [11]. The 'best 
fit' points in their Fig. 15 are 5 m = 60 GeV, m 1/2 = 300 GeV, A = 300 GeV for 
tan/3 = 10 and m = 550 GeV, m 1/2 = 500 GeV, A = 1280 GeV for tan/3 = 50, 
both obtained with m t = 178 GeV and mb(mb) = 4.25 GeV [32]. Both are scenarii 
of stau co-annihilation. The relevant masses, mass differences and the resulting values 
for Q of Ref. [11] (SSARD) are given in Tables 9 and 10 together with the predictions 
from SOFTSUSY 1 . 9, SPHENO 2.2.2 and SUSPECT 2 . 3; we leave out ISAJET 7 . 71 where one 
cannot adjust m^m^). For the point with tan/3 = 10, the Xi an d T\ masses of SSARD, 
which has full 2-loop RGEs, are roughly 2% higher than those of SOFTSUSY 1 . 9 and 
SPHENO 2.2.2. The x^-fimass difference and consequently also Q lie within the values 
of SOFTSUSY 1 . 9 and SPHENO 2.2.2. For the point with tan j3 = 50, however, only SSARD 
has a viable spectrum with a neutralino LSP, while the three public codes get a f\ LSP, 
about 40-60 GeV lighter than the f\ in SSARD. Note also the ~ 10% heavier m A from 
SSARD as compared to the public codes. We can recover a similar AM(x§fi) as [11] with 
SOFTSUSY 1.9 and SPHENO 2.2.2 for A = 1170 GeV. In this case we get m^o = 205- 
206 GeV, m fl = 218-222 GeV, m A = 500-510 GeV and ~ 0.098. However, 1 the fact 
remains that at large tan/3 (and large A ) there are sizeable differences between SSARD 
and the public codes. 

4 Online comparison 

For an easy and user-friendly comparison of SUSY spectrum codes, we have set up a web 
application at 

http : / /cern . ch/kraml/ comparison/ 

Here the user can input mSUGRA parameter points in a web form. The value of the 
top-quark mass is also taken as an input while mb(mb) and a s are fixed to the values hard- 
coded in ISAJET. The mass spectra are then calculated by the latest versions of ISAJET, 
SOFTSUSY, SPHENO and SUSPECT and compared in an output table. The corresponding 
values for Q, Sp, 5a fJ- , B(b — > 57) and B(b — > sp + fi~) are calculated with micrOMEGAs and 
also given in the table. SOFTSUSY is used with the option of full 2-loop running, as in this 
paper. For technical reasons, for the computation of Q a 'static' version of micrOMEGAs 

5 Note that Ref. [11] uses the opposite sign convention for the trilinear A couplings! 
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431.5 


mfj — m^o 


9.7 


8.4 


9.8 


11.7 


m,4 — 2m^o 


181 


193 


196 


196 




0.103 


0.092 


0.109 


0.129 



Table 9: Relevant masses and mass differences (in GeV) and the resulting Q for mo = 
60 GeV, m V2 = 300 GeV, A = 300 GeV, tan/3 = 10, \i > 0, m 4 = 178 GeV and 
m b {m b ) = 4.25 GeV. 





SSARD 


S0FTSUSY 1 . 9 


SPHENO 2.2.2 


SUSPECT 2 . 3 




211 


206 


205 


206 


n 


226 


167 


163 


185 


h° 


117 


116 


116 


116 


m A 


553 


495 


504 


490 


mfj — Tt^o 


15 


-39 


-42 


-21 




0.119 



Table 10: Relevant masses and mass differences (in GeV) for mo = 550 GeV, m\/2 = 
500 GeV, A = 1280 GeV, tan/3 = 50, /i > 0, m t = 178 GeV and m b (m b ) = 4.25 GeV. 

is used which is limited to (co)annihilation channels initiated by Xi2 3> Xii &Ri fi-R, T\, 
and t\. We have checked that this is largely sufficient within mSUGRA. The webpage is 
also useful for comparisons with other spectrum codes and/or programs computing the 
neutralino relic density. 

5 Conclusions 

We have investigated the impact of uncertainties in SUSY spectrum computations on the 
prediction of the neutralino relic density. To this aim we have compared the results of four 
public spectrum codes, ISAJET7.71, S0FTSUSY 1 . 9, SPHENO 2.2.2 and SUSPECT 2.3, in 
the context of mSUGRA. For 'moderate', i.e. not extreme, values of the model parameters, 
we found that the codes in general agree quite well, at the level of few percent, for the 
prediction of the SUSY spectrum. This is also true at large tan (3. 

Nevertheless these small discrepancies can have a large impact on the prediction of 
the relic density of dark matter. We have studied in detail the most important scenarios 
for neutralino (co) annihilation. In the bulk region (although largely excluded by the LEP 
bound on mh), predictions are under control, that is uncertainties are below the experi- 
mental uncertainties of WMAP. In the co-annihilation region, however, the uncertainties 
can easily exceed 30%. Most of this is related to the mass difference between fi and the 
LSP. For this estimate of uncertainties we have used the predictions from ISAJET7.71, 
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SOFTSUSY 1 . 9 and SPHENO 2.2.2. SUSPECT 2 . 3, which only has 1-loop RG running for 
the sfermion mass parameters as opposed to full 2-loop RG running in the other codes, 
typically finds a higher f\ mass. To reduce the uncertainty originating from the spectrum 
calculation to a level below the experimental uncertainty of WMAP, one needs a precision 
in the Xi~fi mass difference of the level of 1 GeV . This corresponds to computing the 
LSP and NLSP masses to per-mille accuracy. Already it has been shown that going from 
2-loop to 3-loop RGE running [33] induces corrections of about that level. 

Similar arguments hold for scenarios where the neutralinos annihilate through pseu- 
doscalar exchange. Typically this means an enhanced coupling to the pseudoscalar, that 
is large tan (5. The critical parameter in this case is the tua — 2m^o mass difference. Al- 
though the codes we compared agree at the level of few percent on the pseudoscalar mass, 
this difference together with the difference in the LSP mass can again add up to 30% or 
more uncertainty in the relic density. To improve the precision of — 2m^o, one not 
only needs to go to higher orders in the RG running but also a more precise treatment 
of the Yukawa couplings, especially of hb, is needed eventually including the full 2-loop 
corrections. At this level also a precise treatement of h T becomes important. Notice, 
however, that the dominant source of uncertainty in is still the present error in the 
extraction of m^rrib). 

Models with non-zero A have usually similar features w.r.t. the relic density as mod- 
els with Aq = 0. The exception is the case of a very large Aq where the t± becomes light 
enough to contribute to co-annihilations. The existing ~ 10% uncertainty in the predic- 
tion of the t\ mass can then lead to order-of-magnitude discrepancies in the prediction 
of the relic density. In particular, ISAJET7.71 predicts a lighter ti than SOFTSUSY 1.9 
and SPHENO 2.2.2, and thus a much lower value for Q. The prediciton of SUSPECT 2 . 3 for 
mi is much below that of the other programs. In fact, in the t\ co-annihilation region 
of SOFTSUSY 1.9 and SPHENO 2.2.2 (but also in the one of ISAJET7.71), SUSPECT 2.3 
does not provide a viable spectrum due to a t\ LSP. This underlines the importance of 
including the full 2-loop RG running in the sfermion masses. 

The picture is however different for extreme scenarios with very large m . These 
are the most difficult models to handle, and large discrepancies in the prediction of the 
spectrum calculators are found. This is especially the case for one of the most important 
parameters for the calculation of the relic density, /i, which determines the masses and 
higgsino fractions of the neutralinos. Predictions for /i can vary by a factor of 2 or more, 
inducing huge order-of-magnitude differences in the relic density. An improvement of 
the situation requires in particular a much more precise computation of the top Yukawa 
coupling. A precise measurement of the top-quark mass, as addressed in [34] , would also 
reduce the uncertainty. Owing to the extreme sensitivity of fi to the exact value of ht 
near the border of REWSB, we consider this region as very unstable. 

We conclude that when using the WMAP bound for constraining mSUGRA models, 
uncertainties from the spectrum computation should be taken into account in addition to 
the experimental uncertainty of Q. For an estimate of the theoretical uncertainties one 
may use the maximal and minimal exclusion curves of different state-of-the-art codes, as 
we have done in this paper. The 5Q obtained this way is comparable to the one obtained 
in [18] by varying the renormalization scale within a given spectrum code. Finally, this 
theoretical uncertainty should also be combined with the uncertainty arising from the SM 
input parameters. 

In parameter regions where SQ originating from spectrum uncertainties is at present 
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larger than the experimental uncertainty from WMAP, more precise calculations are cer- 
tainly desirable to improve the reliability of relic density predictions within GUT-scale 
models. Such improvements will be even more important in view of the precision envisaged 
by the PLANCK experiment. 
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Note added 

After this paper has appeared as a preprint on hep-ph, a new version of SUSPECT, v2 . 3 . 4, 
was published including the 2-loop RGEs for squark and slepton parameters. Owing 
to this improvement, the sfermion masses obtained with SUSPECT 2. 3. 4 agree well with 
those of SOFTSUSY 1 . 9, the two programs now being on the same level in the imple- 
mentation of radiative corrections. In particular, for the t co-annihilation point of Ta- 
ble 8, SUSPECT 2. 3. 4 now gives a viable spectrum similar to that of SOFTSUSY 1.9 or 
SPHENO 2.2.2, with m^o = 143 GeV, m~ tl = 178 GeV and Q = 0.153. This confirms our 
observation of the importance of these 2-loop terms. We note, however, that this does 
not change the 5Q shown in the figures, since SUSPECT 2 . 3 was not taken into account for 
the estimate of uncertainties. 
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